Implementing global Abelian symmetries in projected entangled-pair state algorithms 
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Due to the unfavorable scaling of tensor network methods with the refinement parameter M, new 
approaches are necessary to improve the efficiency of numerical simulations based on such states in 
particular for gapless, strongly entangled systems. In one-dimensional DMRG, the use of Abelian 
symmetries has lead to large computational gain. In higher-dimensional tensor networks, this is 
associated with significant technical efforts and additional approximations. We explain a formalism 
to implement such symmetries in two-dimensional tensor network states and present benchmark 
results that confirm the validity of these approximations in the context of projected entangled-pair 
state algorithms. 

PACS numbers: 75.40.Mg,03.65.Ud 



I. INTRODUCTION 

The density matrix renormalization group (DMRG) 1 
and matrix-product states (MPS)2 have proven to be ex- 
tremely powerful algorithms for one-dimensional quan- 
tum systems. For highcr-dimcnsional systems, however, 
they scale unfavorable with the system size. The reason 
for this is found in the scaling of entanglement entropy, 
which is for many systems governed by the area law. This 
scaling cannot be correctly captured with MPS. 

Other ansatz states have been proposed that by con- 
struction obey the correct scaling of the entanglement 
entropy. Prominent classes of such states are projected 
entangled-pair states (PEPS)^~— and the multi-scale en- 
tanglement renormalization ansatz (MERA)j2s — Just as 
DMRG/MPS, these ansatz states cover the full Hilbert 
space of the quantum problem with a systematic refine- 
ment parameter M. In MPS algorithms, the scaling of 
the computational complexity with this refinement pa- 
rameter is 0(M 3 ), where the number of variational pa- 
rameters grows as 0(M 2 ). In the case of, e.g., PEPS 
on an infinite square lattice, the scaling of computational 
complexity is 0(M 12 ) while the number of variational pa- 
rameters grows as 0(M 4 ). This extremely fast increase 
of computational effort severely limits the attainable M 
to currently about M = 2 ... 8. 

Previous work22 has shown that the accuracy that can 
be obtained with such limited bond dimensions is very 
limited in particular for gapless with a large symmetry 
group. In order to make progress towards controversial 
problems in condensed matter theory, it is therefore nec- 
essary to significantly improve the accuracy of tensor- 
network state methods by reaching larger bond dimen- 
sions M. 

In one-dimensional DMRG calculations, exploiting 
global Abelian symmetries has led to large improvements 
of the accuracy^ Non-Abelian symmetries have also 
been considered.— ~— In the context of two-dimensional 
tensor network state calculations, symmetries have only 
been explored very recently. Parity symmetry (Z2) 



plays a central role in the definition of fermionic tensor 
networks 3 - 8 . - — but has also been shown to be useful for 
spin systems4£ Continuous groups, such as U(l), have 
been used in calculations with the TERG algorithm 4 ^ 
and the MERA. 48 ' 49 A general introduction to the topic 
without numerical results is given in Ref. [50|; Ref. con- 
tains a detailed introduction to U(l) symmetry and its 
use for MERA computations. 

In this paper, we will develop a formalism to imple- 
ment Abelian symmetries into tensor network states. We 
will study the example of infinite projected entangled- 
pair states and numerically confirm the validity of the 
approximation introduced by restricting the structure of 
the tensors. 



A. Projected entangled-pair states 

Let us now turn to a short introduction of projected 
entangled-pair states. We consider a lattice system with 
a tensor-product Hilbert space H = Hi and a product 
basis {\(f>) = \4>i) \(f>2) ■ ■ ■ }■ In order to approximate the 
coefficients c(<j>) of a wave function \^f) = ^ c(<p)\<p) , we 
associate with each site of the physical lattice a tensor of 
rank z + 1, where z is the number of nearest neighbors 
of the site. In this paper, wc will focus on the square 
lattice, where z = 4. A graphical representation is shown 
in Figure [TJ Of these z + 1 indices, one is considered the 
physical index of the tensor with dimension d = dimH;, 
whereas the other ones are auxiliary indices connecting to 
the nearest neighbors with dimension M. The coefficient 
c(4>) is then given as the trace over all auxiliary indices 
in the network. 

To represent a lattice with TV incquivalent sites, usu- 
ally TV different tensors have to be optimized. We can 
however assume that the system is invariant under trans- 
lations by a certain number of sites. Such a state can be 
represented with only few independent tensors and the 
thermodynamic limit can be taken directly. 

PEPS are a highcr-dimcnsional generalization of 



2 




Figure 1. Pictorial representation of a projected entangled- 
pair state (PEPS). Left panel: For the square lattice, a ten- 
sor of rank 5 will be associated with each lattice site. The 
index pointing down connects to the physical system, while 
the other indices connect to neighboring tensors in the state. 
Right panel: The panel shows the PEPS decomposition of a 
coefficient c(</>) for a state |^) = c(<f>i . . . cj)g)\(j)i . . . 4>q) on 
a 3 x 3 square lattice with open boundary conditions. 

matrix-product states. They inherit important proper- 
ties from MPS: i) For M = 1, they are equivalent to static 
mean-field theory, ii) They can capture the entanglement 
properties of gapped systems in the sense that the rank 
of reduced matrices for a block of sites is bounded by 
the exponential of the surface of the block, which al- 
lows the entanglement entropy to diverge with an area 
law. Unlike matrix-product states, however, the exact 
evaluation of expectation values can in general not be 
performed in polynomial time. Therefore, approximate 
methods are required. Several such methods have been 
proposed . 3 i 15 ' 21 i 51 They all have in common that they 
lead to polynomial scaling, yet with large exponents. 

II. SYMMETRY GROUPS 

A. Charge calculus 

In the following, we will be concerned with Hamiltoni- 
ans H with a symmetry group Q , i.e. they commute with 
the elements of some Abelian group Q, [H, q] = \/q G Q. 
This implies that eigenstatcs of H are also eigenstates of 
q. We require: 

• There exists a unitary representation U of the 
group. For q G Q, we have U T (q^ 1 ) = U(q). 

• All representations of the group decompose into a 
direct sum of irreducible representations Vj, which 
can be labelled in correspondence to the eigenvalues 
of some operator g. We will call these labels 

• Consider a state \<f>) G V% and q G Q. We then have 

U{q)\<j>) = v{q,c i )\ct>) (1) 

with v(q, a) G C. 

Examples will be discussed in Section III Bl 

We can classify eigenstates of the Hamiltonian H into 
the irreducible representations of Q. The associated 



labels Ci are then called good quantum numbers (for 
brevity, we will also call them charges or symmetry sec- 
tors). 

As we will see below, only few properties of the irre- 
ducible representations of the group are needed to imple- 
ment symmetric tensor networks. These are intimately 
related to the properties of the eigenvalues defined in 
Eq. ©. 

Charge fusion: Consider two states |o) G Hi, \b) G H2 
with associated quantum number c\ and C2 , respec- 
tively. Their tensor product in "Hi <S> W2 also has 
a well-defined quantum number C3. We thus define 
the fusion of quantum numbers 

ci x c 2 = c 3 . (2) 

This corresponds to a labeling of the tensor product 
of irreducible representations. For the eigenvalues 
this corresponds to 

v(q, cx)v{q, c 2 ) = v(q, c 3 ). (3) 

Identity charge: There exists an identity charge I such 
that c x I = c Vc. This implies v{q,V) = 1. 



Conjugate charge: For each charge c, a conjugate 
charge c exists such that 

c x c = I (4) 

This imples v{q,c) = l/v(q,c). 

We can easily generalize the above to products of 
groups. For Q = Q 1 x (J 2 , the irreducible representations 
are Vy = V} x V 2 , which can be labelled by = (cj, c 2 ). 
These labels correspond to eigenvalues of the operator 
g = (g 1 ® I, I ® g 2 ). The above calculus is then con- 
structed from element-wise operations on the c. 



B. Examples 

An important example is the U(l) symmetry, which 
is present in systems with particle number conservation 
and many spin models. For benchmarking purposes, we 
will apply the symmetric PEPS algorithm to a system 
of spin-i degrees of freedom on the square lattice with 
Heisenberg interaction. This system has an SU(2) spin 
rotation symmetry, which in the thermodynamic limit 
and at zero temperature is spontaneously broken to a 
U(l) symmetry. We will exploit this group and its finite 
subgroups. 
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1. Z 2 

For Hamiltonians that arc invariant under a simulta- 
neous flip of all spins, [ t*} -!•)> ^ e operator 

fe=(-i) Ei<7? =n°? ( s ) 

i 

commutes with the Hamiltonian. A unitary representa- 
tion of the group Z 2 is given by 

U(a) = gl (6) 

with a G {0, f }. Its unitarity follows from the unitarity of 
the Pauli matrix a z . The two irreducible representations 
can be labeled as c = ±. The eigenvalues v(a,c) are: 

i/(0, +) = +1 i/(l,+) = +l (7) 
i/(0, -) = +l i/(l, -) = -l (8) 

The fusion rules therefore are: 

±x± = + ±x=f=- (9) 

This implies + = I and c = c. 

Due to the very simple structure with only two irre- 
ducible representations, the implementation of Z 2 sym- 
metry is particularly easy. 

2. (7(1) 

The most commonly used symmetry in simulations 
with exact diagonalization or DMRG is the U(l) spin 
symmetry, which is given if the operator 

i 

commutes with the Hamiltonian. This is the infinitesimal 
generator of a representation of (7(1), 

U{4>) = cxp (i2n(t>g u{1) ) , (11) 

where <fr G [0, 2n). 

The irreducible representations can be labeled with in- 
teger numbers, c G Z. The i/(0, c) are 

c) = exp {i2-K(f)c) . (12) 

Clearly, 

v(4> : c{)v(4>,C2) = v{<f>,ci+ C2) (13) 

i/(0, 0) = 1 (14) 

u(4>,-c) = vi^c)- 1 (15) 

The charge calculus therefore follows the rules of integer 
addition. The label of the irreducible representations can 
be interpreted as magnetization of the state. Special care 
must be taken when forming the adjoint of a vector or 



operator, since U{4>)YS>) — > (^>\U(—(f>). The Hcrmitian 
transpose of a state in the irreducible representation c 
therefore falls into the irreducible representation c. 



3. l q 

Since for the PEPS, finite groups are easier to deal 
with, we consider finite subgroups of U(l), namely the 
cyclic groups Z q . We define 




which naturally also commutes with the Hamiltonian if 
gu(i) docs. The irreducible representations can be la- 
beled with c G {0, ...,q}, where is the identity. A 
unitary representation is, similar to Z2, given by 

U(a)=gl. (17) 

where a G {0, . . . , q — 1}. The eigenvalues c) are 

v{a,c)=(e i2 f) aC . (18) 

This implies the cyclic property v(a, c+q) = v(a+q, c) = 
v(a,c). The resulting fusion rule is 

Ci x c 2 = (ci + c 2 ) mod q, (19) 

therefore 

c=q-c. (20) 

For taking adjoints, the same consideration as in the case 
of U(l) applies. The implementation of 7L q symmetry 
for q > 2 is more involved than Z 2 since charges are 
not inverse to themselves. The small number of sectors 
however reduces the technical efforts. 



III. SYMMETRIC TENSOR NETWORKS 

A. Definition and contraction of symmetric tensors 

We define a tensor T as a linear map from a tensor 
product of Hilbert spaces to the complex numbers: 

T :Hx®H2® ...®Ur->C. (21) 

Here, R is the rank of the tensor. The elements of the 
tensor are T(v\, V2, ■ ■ •) for Vk G Tik- Equivalcntly, if wc 
choose a fixed basis {b^} in each Tik, we can define a 
tensor as a multidimensional array T^i^..., where the 
indices ik run from 1 to dim%fc and 

m].. !>:..!> ;,...). (22) 
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In this paper, wc are interested in states composed 
of tensors that are invariant under the operations of a 
group. To define this, let q G Q and U k (q) unitary rep- 
resentations in the Hilbert spaces Hk- We then require 

T(U\q)v u U 2 (q)v 2 , ...) = T(«i, v 2 , . . .)■ (23) 

As shown in the Appendix, a tensor element 
T{v\, V2, ■ ■ ■) vanishes unless 



(a) 



Xc k = 

k 



(24) 



where Ck is the label of the irreducible representation 
that Vk belongs to. Colloquially, this can be understood 
as conservation of charge at the tensor. As a direct con- 
sequence, if a fixed basis of eigenvectors of the genera- 
tors is chosen, the multidimensional array T^i^... takes 
a block-sparse form, therefore reducing the number of 
non-zero parameters. 

If we partition the indices to form two groups I± , I 2 , we 
can equivalently express the tensor as a linear operator 



T:Q<)Hk 

keii 



(g) n h 

kei 2 



(25) 



where (®kex 2 Vk)*T(® keIl Vk) = T(v 1 ,v 2 ,v 3 , . . .). We re- 
fer to X\ as in-going and I 2 as out-going indices. In 
a pictorial representation, we will associate arrows with 
the indices. What are the symmetry properties of this 
operator? As shown in the appendix, it commutes with 
the group action. Schur's Lemma then implies that for 
x G (S^fceii m ^ ne inducible representation labeled 
c, Tx is also in the representation c. This is true for all 
possible partitions of the indices. 
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B. 



Tensor contraction 



The steps involved in the contraction of two rank-4 ten- 
sors over two indices are shown in Fig. [2] It is important 
to note at this point that in order to define these opera- 
tions, only the charge calculus introduced in Section [II Al 
is necessary. It is not necessary to know matrix repre- 
sentations of the group in all Hilbert spaces. This will 
allow us to introduce tensor networks with symmetries 
not just on physical, but also auxiliary bonds. 

The steps of the contraction of two symmetric tensors 
are: 

i) Wc first transform the tensors to operators of the 
form (|25J) (Fig. [2J (a)-[2J (b)). The choices of in-going 
and out-going indices are dictated by the indices that 
are being contracted: on one tensor, those indices 
must be the in-going and on the other the out-going 
indices. The resulting operators, which are written 
as a matrix, have a block-diagonal structure. 

ii) The contraction is now equivalent to a matrix mul- 
tiplication. The blocks must be contracted in such 



Figure 2. Pictorial representation of the contraction of two 
rank-4 tensors with 1/(1) symmetry. The steps are explained 
in detail in Section Till Bl 



a way that the resulting tensor still satisfies 
Therefore, we must match blocks in such a way that 



Cin X C ou t — I- 



(26) 



iii) The resulting tensor, Fig. [2] (c), can be converted 
back to the form of Eqn. ([2Tj). 



The conversion between the forms (|2"Tj) and (|2"5| also 
allows the definition of other linear algebra operations, 
such as singular value decomposition, eigenvalue decom- 
position based on the mapping to a matrix. All these 
share the block-diagonal structure. 
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C. Symmetric PEPS 

As a simple example of a tensor network, the construc- 
tion of a matrix-product state invariant under some sym- 
metry group Q is shown in Fig. [3] On each bond i of the 
tensor network, we have a set of charges Cj. For the 
bonds connecting to physical degrees of freedom, this set 
of charges is fixed by the physical Hilbcrt space. In the 
case of a finite MPS with open boundary conditions, the 
set of charges possible on an auxiliary bond is unique and 
has a well-defined physical meaning: if one were to con- 
sider, e.g., a system with particle number conservation, 
the allowed symmetry sectors on each auxiliary bond in 
the construction in Fig. [3] are the possible particle num- 
bers to the left part of the chain. In general, the set of 
allowed charges corresponds to the possible fusion out- 
comes of all physical charges to the left. In a finite sys- 
tem, a quantum number sector can be selected by appro- 
priately fixing the allowed charges at the right end of the 
chain. 

In the case of a PEPS, a unique identification of the 
charges on an auxiliary bond with the fusion outcomes 
of a specific region cannot be made. It is therefore not 
possible to determine uniquely which symmetry sectors 
must be kept on the auxiliary bonds. While for finite 
groups, it is usually computationally possible to allow 
all charge sectors, some choice has to be made in the 
case of infinite groups. It will therefore be one of the 
main purposes of this paper to verify that i) for finite 
and infinite groups, one obtains a good approximation 
to the ground state by using a PEPS constructed from 
symmetric tensors, ii) for infinite groups, a reasonable 
approximation is obtained for computationally feasible 
choices of the symmetry sectors. 

To understand the nature of the approximation intro- 
duced by truncating the set of allowed quantum numbers, 
consider the expansion of a state |\&) = J2\<p) C (0)I0)- Us- 
ing a tensor network ansatz amounts to representing all 
coefficients c(4>) by a trace over a tensor network, which 
will represent the low-entanglement subspace of the full 
Hilbcrt space efficiently In principle, all basis states are 
allowed and could have non- vanishing weight. Imposing 
restrictions on the quantum numbers, on the other hand, 
amounts to a restriction on the allowed basis states: the 
sum does not run over the full basis {\(p}}, but only a 
subset of states compatible with the allowed quantum 
numbers. 

In addition to the charge sectors on each bond, the 
number of states in each sector has to be chosen. In 
principle, this could differ between all sectors on one 
bond and between bonds. The situation becomes more 
involved since even for a translational invariant PEPS, 
several different sets of charges have to be considered, 
as shown in Fig. 01 For the purpose of this paper, we 
make the simplification that we choose the charges to be 
the same on all equivalent bonds of the lattice and the 
environment states. Additionally, for the case of finite 
groups, we choose the number of states in each sector 




c 
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Figure 3. End of a matrix product state invariant under some 
symmetry group Q. \<f>i) denote physical states in the local 
Hilbert space %ioc- By C, we denote the set of charges asso- 
ciated with sectors in Hi oc , and by C n the set of charges as- 
sociated with sectors in (^)™ =1 Hioc- The first auxiliary bond 
to the left simply carries the physical charges of the first site. 
For the second bond, all possible fusion outcomes of charges 
on the first auxiliary bond with the physical charges have to 
be considered. This can be continued up to the middle of 
the chain, such that each auxiliary bond carries the possible 
combinations of charges to the left. Joining such a state with 
its reflection will yield a finite symmetric MPS. 

the same for all equivalent bonds. 

Two points require special attention when applying 
symmetric PEPS to infinite lattices: i) On infinite lat- 
tices, the ansatz is restricted to states that globally fall 
into the sector of the identity charge. For example, using 
the U(l) symmetry of a spin-i system described in Sec- 
tion lll B~2) only states with vanishing total magnetization 
can be studied. In the case of particle number conserva- 
tion, an appropriate choice of charges would have to be 
taken to enforce the desired filling fraction. On finite lat- 
tices, however, selecting specific quantum number sectors 
is possible also in the PEPS construction by adding an ex- 
ternal bond carrying the total charge of the system to one 
of the tensors that make up the PEPS, ii) Since our con- 
struction assumes that the state has a well-defined global 
quantum number, systems that spontaneously break the 
symmetry that is being exploited cannot be represented. 
If the possibility of a spontaneous symmetry breaking is 
present in the system being studied, the results should 
therefore be checked against calculations without enforc- 
ing the symmetry. 

D. Implementation 

In this section, we will outline a few details of our im- 
plementation. The most important operations on tensors 
include i) contraction, ii) singular-value decomposition, 
and iii) eigenvalue decomposition of tensors. In particu- 
lar for the last two operations, very efficient implemen- 
tations exist for matrices and it is advisable to make use 
of these. The contraction could in principle be implc- 
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Figure 4. Corner of a symmetric PEPS state with an environ- 
ment as in the directional corner transfer matrix method.— 
Here, the blue circles denote tensors T; of the ansatz state and 
their conjugate T* , the red squares denote tensors of the cor- 
ner transfer matrix and the orange circles represent single-site 
operators acting on the physical index of the PEPS tensor. In 
the infinite case, in general there are three sets of charges in- 
volved: i) the physical charges C p h ys carried on the black, solid 
lines in the figure, ii) the auxiliary charges C aux on the blue, 
dashed bonds, and iii) the charges carried on the bonds of the 
environment Ccomcr- This reflects the three independent bond 
dimensions involved in a PEPS: the physical dimension d, the 
bond dimension M and the environment dimension \. Usu- 
ally, M > d and \ ~ M 2 ■ Therefore, C p h ys C C aux C C CO mer- 
In principle, all charges could depend on the location in the 
PEPS or the environment. 



merited directly as a summation, however it turns out to 
be favorable to map it to matrix multiplication and make 
use of existing, optimized implementations. 

The most important operation therefore is the map- 
ping between a symmetric tensor and a block-sparse ma- 
trix, i.e. between the form ([2"Tj) and (|2~5j) . with sev- 
eral tensor indices grouped to form the left and right 
indices of the matrix. Since this operation, like a ma- 
trix transpose, scales like O(N), with N the number 
of elements in the matrix, it is subleading compared to 
contraction, singular- value decomposition, etc., which all 
scale roughly as (D(N 3 / 2 ). 

For the conversion between tensors and matrices, one 
could either calculate the correspondence between the lo- 
cation of an element in the tensor and in the matrix on 
the fly or compute it once and store it in memory along 
with the tensor (precomputation). This is explained in 
some detail in Ref. H^. We choose not to use precom- 
putation for several reasons: a) the overhead in memory 
usage may be significant, b) memory bandwidth is one of 
the bottlenecks of tensor-network state simulations and 
should therefore be minimized, c) the structure of ten- 
sors, in particular for groups such as U(l) where the 
number of sectors is chosen dynamically, may vary be- 
tween each iteration of the algorithm, d) the overhead of 
calculating the tensor structure on the fly is negligible if 
implemented efficiently in a compiled language such as 
C++. 
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Figure 5. The first panel shows how the relative error in 
the energy decreases as M is increased, for different choices 
of the symmetry group. The second panel shows the same 
data versus the number of variational parameters in the state 
(note the logarithmic scale on both axes). Clearly, the fact 
that the relative errors are similar between symmetry groups 
for a given M shows that reducing the number of parameters 
and the computation time by using larger (finite) symmetry 
groups does not lead to any loss in accuracy. 
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Figure 6. Relative error in the ground state energy of the 2d 
Heisenberg model with a [/(l)-symmetric PEPS, as a function 
of i) the total bond dimension, ii) the number of parameters. 
The choice of sectors and dimensions is shown in Tabled] We 
show results with Z2 symmetry for comparison. 



IV. RESULTS 

It has been demonstrated in Ref. that the spin-i 
Heisenberg model, 



(27) 
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where Si is the spin-^ operator at site i and the summa- 
tion runs over pairs of nearest neighbors, is a difficult test 
case for tensor network methods. As in Ref. HH, we will 
work on an infinite square lattice. This is due to strong 
fluctuations around the Neel state, which reduce local 
magnetic moments significantly. We will use the Heisen- 
berg model as a benchmark case here and compare to 
precision Monte Carlo calculations^^ 

All results in this section were obtained using the sim- 
plified update described in Refs. [ID and |H. In this up- 
date scheme, long-range correlations are effectively taken 
into account by introducing weights on the auxiliary 
bonds of the PEPS. Imaginary time evolution is then 
performed locally, determining new tensors and weights 
bond after bond. While no formal justification can be 
given for the weights, the accuracy of the algorithm ap- 
plied to systems away from criticality turns out to be 
only slightly less than an update scheme that takes cor- 
relations into account more rigorously. The advantage of 
the simplified algorithm lies in the much better perfor- 
mance and robustness against numerical instabilities. 

For the imaginary time evolution, a Suzuki- Trotter de- 
composition of the evolution operator has to be per- 
formed. Since numerical errors do not accumulate in 
imaginary time evolution, we can reduce the time dis- 
cretization during our simulation and completely sup- 
press discretization errors. To extract expectation values, 
we use the directional corner transfer matrix approach of 
Ref. lU We use an ansatz with 4 independent tensors in 
a 2 x 2 unit cell. 



A. Finite groups 

The results we obtain for the Heisenberg model with fi- 
nite symmetry groups Z2 and Z3 are shown in Fig. [5] as a 
function of the total bond dimension on the bonds of the 
state and as a function of the total number of variational 
parameters of the state (note the logarithmic scale in this 
case). For comparison, we show results obtained with a 
non-symmetric PEPS, but with the same simplified up- 
date scheme. We choose the number of states equal in 
each sector, hence M = q ■ n. We also make the same 
choice on all bonds of the PEPS. We keep up to 36 states 
in the renormalization of the corner transfer matrix. 

For n > 1, that is with a non-trivial dimension in 
each symmetry sector on the auxiliary bonds, the ener- 
gies obtained with the symmetric PEPS are comparable 
to those obtained without symmetry for the same bond 
dimension. This demonstrates that the approximation 
introduced by restricting the structure of the tensors is 
valid and does not affect the accuracy. Since all ma- 
trix operations decompose into q blocks, we can expect a 
speedup of 0(q 3 ) of the algorithm. In terms of the num- 
ber of variational parameters, a significant improvement 
is achieved: with Z2 symmetry, only half the number of 
variational parameters is necessary. With Z g symmetry, 
the reduction is even stronger. This may be advantageous 
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Comparison 
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6 


2-2-2 


2048 


5184 
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1-2-2-2-1 


4800 


16384 


5 
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1-2-3-2-1 
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2-2-2-2-2 


10240 
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5 


11 


2-2-3-2-2 


15680 


58564 


7 


14 


2-2-2-2-2-2-2 


28672 


153664 



Table I. The table shows the choices of symmetry-sector di- 
mensions for the U(l) symmetry in the Heisenberg model. 
The columns contain i) the number of sectors associated with 
quantum numbers S z = — (n — l)/2. . . (n — l)/2, ii) the to- 
tal bond dimension, iii) the size of each sector, iv) the total 
number of parameters of the state, v) the total number of 
parameters in a Z2-symmetric state with the same total bond 
dimension. 



particularly if a direct energy minimization algorithm is 
applied instead of the imaginary time evolution. 

In some cases, the energy of the symmetric state falls 
below the energy of the non-symmetric state. This must 
be attributed to trapping in local minima, which seems 
more likely in the case of a non-symmetric PEPS with 
more variational parameters. 



B. 17(1) 

While for the finite groups considered so far, we could 
simply keep all allowed sectors of the symmetry on each 
auxiliary bond, some choice must be made for the in- 
finite group U(l). Furthermore, we have to choose the 
dimension within each symmetry sector - due to the large 
number of sectors, it is generally not efficient to keep it 
the same in all sectors, as we did for finite groups. How- 
ever, given the fast growth of computational cost with 
the bond dimension, only a few choices are possible. The 
choices we considered are listed in the table in Fig. |TJ It 
should be noted that for equal total bond dimension M, 
a state with more symmetry sectors of smaller dimension 
is computationally less expensive since all matrix com- 
putations can be split into more blocks. This allows us 
to study states with very large bond dimension up to 
M = 14, which would be intractable otherwise. 

Results obtained with the above choices are shown in 
Fig- El The accuracy for a given bond dimension is worse 
than with the finite group Z2; even for the very large 
bond dimensions studied with U(l) symmetry, the ac- 
curacy does not reach the level of the finite symmetry 
groups. This is a clear signature that the approximation 
we made by imposing a J7(l)-symmetric structure on the 
tensors and picking only a few allowed sectors of the sym- 
metry limits the accuracy of the simulations. One has to 
keep in mind, however, that the number of variational pa- 
rameters is reduced much more strongly than in the case 
of finite groups, as shown in the last column of Table [IJ 
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V. CONCLUSION 

Wc have explained a formalism for introducing Abclian 
symmetries into tensor network state algorithms. The 
formalism relies only on fusion properties of irreducible 
representations and is therefore easily applied to a large 
class of symmetry groups. The formalism can be applied 
to any tensor network state algorithm; for this paper, we 
have restricted ourselves to simulations with projected 
entangled-pair states. 

Since the implementation requires additional approxi- 
mations, benchmark calculations confirming the validity 
of the approach are required. This is particularly im- 
portant in the case of U(l), where restrictions on the 
allowed quantum numbers have to be introduced. In or- 
der to assess the validity, we have applied our method to 
the spin-i Heisenberg model on an infinite square lattice 
with the symmetry groups Z 9 for q = 2, 3 and U(l). 

Our results for the finite groups show that no accuracy 
is lost due to the symmetric decomposition of tensors. 
At the same time, the number of variational parameters 
and the computational effort is significantly reduced. We 
therefore expect that exploiting these symmetries will be- 
come very useful in the context of tensor networks states. 
As a future direction of research, a decomposition where 
the number of states in each symmetry sector is not equal 
could be considered, which may lead to even better ac- 
curacy for a given bond dimension. 

In the case of the continuous symmetry group U(l), 
we were able to achieve much larger bond dimensions. 
Nevertheless, the accuracy does not reach the level that 
can be obtained with finite symmetry groups. We expect, 
however, that if a sufficiently large number of symmetry 
sectors is taken into account, the accuracy will eventually 
become comparable to the non-symmetric case. Further 
research is required to understand this, in particular how 
this behaves for different models such as bosonic models 
with particle number conservation. Also, a scheme that 
automatically picks the relevant symmetry sectors on the 
auxiliary bonds and in the environment tensors without 
strong dependence on the initial state may be very useful. 

We acknowledge useful discussions with G. Vidal. Sim- 
ulations were performed on the Brutus cluster at ETH 
Zurich. 



VI. APPENDIX 

In this Appendix, we show in detail some calculations 
relevant to the discussion in Sect. lIII "AI in particular sym- 



metry properties of a tensor of the form (12 1 j) . 



To simplify the notation, wc consider the case of a 
tensor of rank 2, where the Hilbert spaces are taken to 
be equal: 

T:H®U^<£ (28) 

f-.n^n, _ (29) 

where for x,y € H we have T(x, y) = y^Tx. Also, let U 
be a unitary representation of the group in H. 

We would like to show the following equivalence: 

(a) [T, U] = (30) 
(6) &T(Ux,Uy) = T(x,y). 

First, we show that (b) follows from (o): 

T{Ux,Uy) ^y^U^TUx (31) 
= y^U^UTx 
= T(x,y) 

Secondly, we show that (a) follows from (6): 

y f [T,U]x = y^TUx-y^UTx (32) 
= T(Ux,y)-T{x,U*y) 
= T(U'<Ux,U" l y)-T{x 7 U'<y) 
= 



The above can easily be generalized for all operators 
of the form ([231) . 



Wc now want to show the validity of (flH)) . Let q 6 Q. 
Then, 

T(U(q)v u U(q)v 2 ,...) (33) 
= T(i/(q,c 1 )vi,i/(q,c 2 )v 2 , ■ ■ ■) 
= i/(q, ci)v(q, c 2 ) ... T(vi,v 2 , ■ ■ •) 
= "(<l, X c k )T(vx,v 2 , . . .) 
= T(vi,V2, ■ ■ ■) 

where we have used Eqns. ([23]) and ([3]). Therefore, 
v(q, X Cfe) = 1 and X = I. 
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